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This paper is concerned with the design and analysis of a fuhy adaptive eigen- 
. value solver for linear symmetric operators. After transforming the original problem 

■ into an equivalent one formulated on ^2, the space of square summable sequences, 

• the problem becomes sufficiently well conditioned so that a gradient type iteration 

can be shown to reduce the error by some fixed factor per step. It then remains to 
realize these (ideal) iterations within suitable dynamically updated error tolerances. 
It is shown under which circumstances the adaptive scheme exhibits in some sense 
^ ' asymptotically optimal complexity. 
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1 Introduction 



O ■ 1«1 Background 



d d 



^ I In a Gelfand triple Ti. ^ X ^ Ti.' of Hilbert spaces, where Ti, is densely embedded in X, 
and the dual pairing induced by the inner product (., .) of X, let 



c-.n^n' (1.1) 

be a linear operator that takes Ti. onto its normed dual Ti'. We wish to find, under certain 
assumptions on £, an eigenpair (A, n) of the problem 

Cu = \8u (1.2) 

corresponding to the smallest eigenvalue A. Here S : u ^ {■,u) canonically embeds Ti in 
H'. Of course, (11.21) is to be understood in the weak sense, i.e. 

{v, Cu) = X{v,Su), veH, (1.3) 



*This research was supported in part by the Leibniz Programme of the DFG, by the SFB 401 funded 
by DFG, the DFG Priority Program SPP1145 and by the EU NEST project BigDFT. 
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In view of the above choice of the dual pairing (■, •) we shall from now on simply write 
(w, w) instead of {v, £w) for v,w eH. 

For typical applications, one can think of X as an L2-space over some domain Q and 
7i as a Sobolev space iJ* of positive real order t (or as a closed subspace of a Sobolev 
space, determined e.g. by homogeneous boundary conditions). 

We shall be concerned with the case that the operator C is symmetric, i.e. 

{Cv,w) = {v,Cw), for all v,w (1.4) 
and bounded and strongly positive, i.e. there exist positive constants cl,Cl such that 

CL\\v\\n < {Cv,v) =: \\v\\l < Cl\\v\\1^, veH. (1.5) 

Note that if the (real) spectrum of C is bounded from below by a negative number, C 
can be shifted by a suitable parameter /i such that C — fiS satisfies (II. 5p . With this slight 
modification, this assuption is met by e.g. the electronic Schrodinger operator [Re Si] . one 
particle Schrodinger operators in M.^ with certain potentials [Weidl IHiSij and eigenvalue 
problems for strongly elliptic (satisfying a Garding inequality) differential operators on 
bounded domains as well as by eigenvalue problems for Fredholm integral operators. 

Note that (11. 5p also implies that £ is a norm-isomorphism from 7i onto H.' whose 
condition is bounded by Cl/cl. 

We shall exclusively deal with eigenvalue problems for which the infimum of the spec- 
trum is an isolated eigenvalue A. For simplicity we shall actually assume that 

. {Cu,u) 
A := mr — — 

is a simple eigenvalue with corresponding eigenvector u, and that that the rest of the 
spectrum is bounded from below by A > A, which means that 



{v,v) 



> A holds for all veH with {Cu,v) = 0. (1.6) 



The conventional approach is to discretize (11. 2p . e.g. by finite elements or finite differences, 
which gives rise to a finite dimensional discrete problem 

AhUh = fiChUh (1.7) 

where, for a given basis of the trial space. Ah, Ch are the corresponding stiffness matrix 
of C and the mass matrix, respectively, and u/j denotes the coefficient vector of the ap- 
proximate eigenvector. Now the issue becomes to solve (I1.7P efficiently within a suitable 
accuracy tolerance associated with the discretization. There is a vast literature on this 
issue, see e.g. Parlett |Parl] . Golub/Van Loan [GvLj . which are standard text books in 
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numerical linear algebra, or Chatelin [Chat], Babuska-Osborn |BabOj concerning Galerkin 
discretizations and Dyaknov |Dyak| , Knyazev/Neymair [BPK] for preconditioned iteration 
schemes. See also |RSZ] for further references and comments. 

In this paper we follow a different line. While the subsequent analysis will apply to 
finite dimensional problems of the form (II. 7p as well, our primary interest is to avoid the 
separation of the discretization and solution process. Rather, we will design an abstract 
iteration scheme that solves the original infinite dimensional problem within a given ac- 
curacy tolerance, where the iterates of this scheme are taken, in principle, from all of Ti. 
We will then show how to realize these updates (up to a perturbation) in a finite dimen- 
sional setting, particularly directing our attention to carrying out this task at possibly 
low computational cost. After completion of this work, we became aware of related work 
in |GG] where an adaptive finite element scheme for elliptic eigenvalue problems is shown 
to converge without giving complexity estimates though. 

In principle, such infinite dimensional iterations can be formulated in different ways, 
focussing either on the convergence of Rayleigh quotients or of eigendirections. While the 
first option - although with the same motivation as in the present work - has been adopted 
in |RSZj . we address here the algorithmic realization of the second option. In fact, in 
|RSZ] , from a somewhat different perspective, we focus on convergence of preconditioned 
iteration schemes per se, whereas in the present paper we develop and analyze a convergent 
adaptive algorithm. Our analysis given here provides (a posteriori) criteria for adaptive 
updates and yields complexity estimates that prove the (asymptotic) optimality of the 
scheme. 

We note that problems of the kind (11.31) may be treated by means of inverse iteration 
or better by Davidson or Jacobi Davidson type methods as well. These approaches require 
the solution of a linear system in each iteration step. For this purpose one may apply the 
adaptive solution strategies proposed in |CDDll ICDD2t ICDD3] based again on various 
types of (infinite dimensional) iteration schemes. In contrast to this strategy, in the present 
approach we will intertwine the outer loop, e.g. inverse iteration or Jacobi-Davidson and 
the iterative solver in the inner loop by updating the Rayleigh quotient in each inner 
iteration step. In fact, the present approach can be viewed as a preconditioned steepest 
descent method for minimizing the Rayleigh quotient. At any rate, one should stress that 
for large systems of linear equations the use of iterative methods would be inevitable, 
anyway. 

As in Jacobi Davidson or in Krylov space methods, one can also apply subspace 
acceleration techniques to improve the present scheme by computing Ritz values of a rather 
small system. This requires the computation of scalar products like those appearing in the 
Rayleigh quotient. Since we deal with infinite matrices in the present paper we can only 
compute these values approximatively, see also Section 14.31 for further details. Regarding 
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these possible improvements, we content ourselves here with brief indications and further 
studies are required. 

We proceed now with showing how the problem (11.21) is transformed into an equivalent 
one formulated over a sequence space which is, however, in some sense better conditioned. 
To this end, we begin with collecting a few facts that are relevant for carrying out this 
program. Moreover, this will motivate the assumptions made in subsequent sections. 

1.2 Transformed Problems 

Instead of projecting the problem (11.21) to a (fixed) finite dimensional subspace Tih C Ti 
of 7i , we shall transform (11.21) into an equivalent problem defined on the sequence space 
£2(^) of square summable sequences of some possibly infinite index set X endowed with 
the inner product 

(v, w) = ViWi, for all v, w G ^2(^) 

iGX 

and induced norm || ■ ||. Note that we use the same symbols for the inner product on £2(2^) 
as for the dual pairing on Ti! x Ti. 

The key ingredient is a suitable Riesz basis \1/ = {ipi : i G X} of 7Y, i.e. there exist 
positive constants Cq,^C^, such that 

c^,\M<\\^Vi%IJi\\n<C^\\w\l for all v G £2(2:). (1.8) 
Then (II. 2p is equivalent to 

Au = ACu, (1.9) 

where 

A := ^,))^^^.^^ , C := ((V^,, ^.))^_^.^^ . (1.10) 

For the operator A the properties (11.51) of C together with (11.81) imply (see e.g. jP] ) 
that 

||A||,,^,, < ClCI, \\A-%,^,, < cl\f, (1.11) 
for the problem on -^2(2^), which in turn means that 

cl4i|v||2< l|v||i<CLC|||v||^ (1.12) 

where || ■ ||^ := (A-, ■). Thus, A-ellipticity is equivalent to bounded invertibility of A on 
i^iX)- One could say that the transformation has a built-in preconditioning effect: The 
original L has been transformed to an operator A which is well conditioned on ^2(2^) 
according to (11.121) . a fact that will play a crucial role in solving (11.91) numerically. 

As for the properties of C, in the above setting the operator C stems from the inner 
product on X and is therefore symmetric, bounded and positive definite. The boundedness 
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of C as a linear operator on i'2(X) results from the continuous embedding of 7i in Af. 
However C is typically not coercive on £2(2^)- 

The role of C is further illuminated when specifying the setting slightly to the following 
situation which may actually serve as a guiding example. In fact, consider again the 
example Ti = H*, X = L2. Here, a suitable choice of the basis \1/ is a wavelet-type basis. 
In fact, such wavelet bases are available by now for a wide range of domains and exhibit 
an important property, namely that a scaled version of \1/ is also a Riesz basis for the 
pivot space X . That means, setting 

D := diag(rfi := {^u^i)^'^ = : i ^ X), 

the collection 

*° := {rf-Vi : ^ e X} 

satisfies 

4l|Dv|| < ||5^i;,^i||;,<C^||Dv||, Dvg£2(X), (1-13) 

where d^,C'^ are again fixed positive constants. It is well-known that for X = L2 and 
T-C = one has di = ~ 2^*1*', where |i| denotes the dyadic level of the 

wavelet ipi. Moreover, the matrices A and D~^CD~^ are spectrally equivalent, i.e. there 
are constants c*,C* such that 

c*||v||d-icd-i < ||v||a < C*||v||d-icd-i, vg£2(X). (1.14) 

Remark 1. In both of the above settings, C is symmetric positive definite and bounded 
but typically not coercive on £2(1) ■ However, C is coercive on the space {x : Dx G ^2(X)} 
equipped with the proper norm in the second example. 

Note also that if the basis functions ipi are pairwise orthogonal with respect to the pivot 
inner product {■,■), then C is actually a diagonal matrix. 

The spectrum of the generalized eigenvalue problem Au = ACu coincides with the 
original spectrum of ( 11. 2p . In particular, since C is positive definite (but not necessarily 
coercive) on £2(1), 12.11 is equivalent to C^^/^AC^^^^y = Ay. The minimal eigenvalue is 
then given by 

• • (C-^/^AC-VV,y) . (Ax,x) 

< A = mm— — = mm ; = mm— -, (1-15) 

V {v,v} y (y,y) x (Cx,x) 

so that the minimal eigenvalue of the transformed problem equals the one of the original 
problem and A is simple. Denoting by u G ^2(X) the corresponding eigenvector, one has 

(Av,v) 



(Cv,v) 



> A, for all ve£2{I) with (Av, u) = 0. (1.16) 



The properties collected above will guide us later when formulating the precise condi- 
tions all subsequent developments will be based upon. 
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1.3 Objectives and Layout 



The above considerations indicate how to arrive at transformed equivalent eigenproblems 
on sequence spaces (11 .Qp that are better conditioned in the sense that the spectrum of 
the resulting matrix A is enclosed in a bounded interval of the positive real semi-axis. In 
particular, the resulting ellipticity on certain orthogonal complements will be seen to give 
rise to iteration schemes that reduce the error of the approximate smallest eigenvalue and 
also the deviation of the corresponding approximate eigenspace from the exact one by 
at least some fixed factor less than one. This is the case for the full infinite dimensional 
problem as well as for any finite dimensional discretization resulting from projecting onto 
spans of subsets of the Riesz basis 

Once the fixed error reduction has been established for the full infinite dimensional 
problem, the principal idea is to mimic this ideal iteration on the infinite dimensional 
problem (that still contains all information) by ultimately carrying out these iterations 
only approximately. One then faces the following two tasks: 

(i) Find appropriate stage dependent tolerances within which each iteration is to be 
solved so as to still guarantee convergence to the exact solution; 

(ii) Devise numerical schemes that realize the approximate application in the involved 
operators within the given target accuracy at possibly low computational cost. 

In this paper we explore the potential of such an adaptive solution strategy on a 
primarily theoretical level. Our central objective is to analyze the intrinsic complexity 
of accuracy oriented eigenvalue computations. In particular, we shall show under which 
circumstances such a scheme has asymptotically optimal complexity in the following sense: 
If the solution u belongs to some approximation space A^, which means that terms in 
the expansion of u suffice to approximate m in 7i within accuracy of order A^~'^, then the 
computational work for realizing a target e remains bounded by a fixed multiple of e"^/*, 
when e tends to zero. 

While the present work is certainly inspired by prior work on adaptive solution schemes 
for operator equations (see e.g. |CDD2t ICDD3] ). there are some noteworthy differences. 
On one hand the iteration is nonlinear, on the other hand, one approximates coefficient 
arrays that are only determined up to normalization. 

The layout of the paper is as follows. After summarizing the fundamental properties 
of the problem in Section [2], we formulate in Section [3] an iterative scheme for the compu- 
tation of the smallest eigenvalue of the infinite dimensional problem (II. 2p (or equivalently 
(II. 9p ) and a corresponding eigenspace. Moreover, we analyze the convergence of this ide- 
alized iteration making essential use of the ellipticity of A — AC on certain orthogonal 
complements. Section H] is devoted to the numerical realization of the idealized scheme 
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based on the approximate realization of residuals. We first show under which circum- 
stances and for which accuracy tolerances such perturbed schemes converge to the exact 
solution of (11.91) . Then, in a second step, we analyze the complexity of such schemes un- 
der the assumption that the involved operators A, C are in a certain sense quasi-sparse. 
We note that this property is known to hold for a wide range of operators. The main 
result is that the proposed adaptive scheme exhibits in some sense asymptotically opti- 
mal complexity. While the underlying adaptive scheme is to be viewed as one possible 
prototype realization with certain asymptotic properties, we conclude the section with 
indicating ways of quantitative improvements. Finally, in Section [5] we briefly indicate a 
typical application background. 



2 Problem formulation 

Motivated by the above considerations, we will in the remainder of this paper restrict our 
treatment to determining the smallest eigenvalue A and the associated eigenvector u of 
the generalized eigenvalue problem 

Au = ACu, (2.1) 

formulated on the infinite dimensional space ^2(^)- As in the previous examples we will 
only consider the case where A, C are symmetric, positive definite matrices defined on 
£2(2^) endowed as above with the norm ||x|p := Xliex l-*"*!^ ~- (^'^)- Let us again stress 
that the index set X could (and actually will) be infinite. 

Further, we will impose the following conditions on our problem (12. ip reflecting the 
properties we have identifled in the framework discussed before. 

Property 1. There exist positive constants 7, F such that 

7 ||x|p < ||x||i < F ||x|p, (2.2) 
where \\ ■ ||a := (A-,.-)2 denotes the K-( energy) -norm, see 

Property 2. For the minimal generalized simple eigenvalue A of \2.1[ there holds 

(Ax, x) 

» < i = "rW^- ''^-'^ 

while there exists a A > X one has for all x with ( Au, x) = 

(Ax, x) . , , , 
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Remark 2. In the framework discussed in the previous section it was already observed 
that C is bounded. Note that this actually follows from Properties 1, 2, namely one has 
forCc := r/A that 

l|x||c < Cc ||x||. (2.5) 

On the other hand, it follows from ^2.2^) together with U.15\) that C is coercive on all 
of C-2 if o-nd only if C is bounded with respect to the X -inner product, a condition that is 
usually not fulfilled, for instance when C is a differential operator. 

3 Basic Algorithm - A Richardson- style method for 
calculating the smallest eigenvalue 

In this section, we shall now formulate an "ideal iteration scheme" for (12.11) . This scheme 
and the tools used to prove convergence will be utilized later in the analysis for our 
adaptive algorithm to be introduced in Section |H therefore, the convergence analysis will 
be given in detail. Let A and A be lower and upper bounds for A and A, respectively. Then 
the iteration reads as follows. 



3.1 Basic Algorithm: 
MINIT 



Require: initial value xo G ^2(2^), ||xo|| = 1, 

a = 2((l-A/A)7 + r)-i 

Iteration: 

For n = 0, 1, . . . do 

Calculate the Rayleigh quotient A*^") = ^^^"'x") 

Calculate the residual r„ = Ax„ — A'-"-'Cx„. 

Let x„+i = x„ - ar„; 

Normalize x„+i = (x„+i, x„+i)"22„_|_j. 
endfor 



This algorithm resembles a preconditioned Richardson iteration x^+i = $("^x„ with 
a stage dependent iteration matrix 

:=I-a(A- A(")C) 

depending on the ra-th iterate x„ and a relaxation parameter a. On the other hand, 
defining the (generalized) Rayleigh quotient A(x) = one has VA(x„) = (cx^ x ) 

so that a steepest descent step with stepsize /3„ takes the form 

\'--^X,j, X^y 
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Thus, the above iteration could be viewed as a steepest descent method for the general- 
ized Rayleigh quotient. This iteration scheme is also known as a preconditioned inverse 
iteration scheme PINVIT, see e.g. |KNt IRSZ] . Note that the preconditioning is already 
provided through the formulation as a well-conditioned problem in £2(2^)- 

This algorithm can be improved by subspace acceleration techniques. For instance, 
the Rayleigh quotient might be minimized according to 

x„+i := argmin {A(v) : v G span {x„, r„}, | |v| | = 1}. (3.1) 

This corresponds to optimal line search in a steepest descent step. However, in order 
to have a technically simpler exposition of the principal mechanisms of our approach 
we confine the subsequent discussion first to the above simple Richardson type iteration 
and will take up possible improvements along with computational consequences later in 
Section 14. 3[ 

In the remainder of this section, we will analyze the properties of Algorithm MINIT 
and in particular prove convergence to the smallest eigenvector, provided the starting 
vector lies sufficiently close to the target u. The main results will then be compiled in 
Theorem [T] at the end of this section. 

3.2 Error reduction 

We shall show that for suitable damping parameters a the above (ideal) iteration (on the 
full infinite dimensional space -^2(2^)) provides approximations to the searched eigenpair 
with a guaranteed fixed error reduction per step. 

3.2.1 Preliminary considerations 

Setting 6n '■= 5(x„) := x„, — u for the normalized iterates x„, we will show that the error 
components that are orthogonal to u tend to zero in the Euclidean i2{T)-noTm. Denoting 
by Px = l^^u the £2-orthogonal projection of x into the eigenspace Uq containing the 
eigenvector u, the orthogonal error component is given by 5^ = (I — P)5„, we introduce 

||5^|| = ||(I-P)5„|| = ||(I-P)x„|| =: sinZ(x„,u). (3.2) 

Note that does not depend on the normalization of u. We shall now show that we 

have an at least linear reduction of the orthogonal error component of the iterates Xn of 
dSlD, i.e. 

llClll <ell^nll, (3.3) 

with ^ < 1, provided the initial vector has a sufficiently small angle with the solution u 
corresponding to the smallest eigenvalue A. 
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Towards this end, we first note that x„± r„, so that we have for the non-normalized 
iterates 

||x.„+i|p = ||x„ - ar„|p = ||x„|p + a^||r„|p > ||x.„|p = 1. (3.4) 
Hence, the error of the £2-iiormahzed iterates x„+i can be estimated by 

ll^^+ill = ^^||(I-PR+i|| < ||(I-P)x„+i|| =: \\SI^,\\. (3.5) 

I l^n+l| I 

To go further, we decompose as follows 

ll^n+ill = ||(I-P)$(")x„|| = ||(I-P)((I-a(A-A(")C))xJ| 

< 11(1 -P) (I-a(A-AC))5„ || + ||a(A(")-A)(I-P)Cx„,||, (3.6) 

where we have used that 

Mu := (A - AC)u = . 

Moreover, using the fact that rangeP C ker M so that 

M = M(I-P), (3.7) 
we conclude, upon using (I — P)^ = I — P and 

(A - XC))Sn = MSn = M(I - P)dn = Md^ , 

that 

< II (I-P) (I-aM)5;[ ll + Cc a |A(") -A|. (3.8) 

We now estimate the two parts in the right hand side of (13. 8p separately in the next two 
subsections. 

3.2.2 Linear part of the iteration scheme 

The main goal of this subsection will be to show the following property of the iteration 
matrix, from which in Lemma HI an estimate for the left part of (13. Sp will easily be derived. 

Lemma 1. Under the assumptions stated above, the matrix 

M:=A-AC (3.9) 

is bounded and i2-elliptic on the set 

\/o = {xG£2(X):(x,u) = 0}, (3.10) 
that is, there exist 6,Q > such that 

^ ||x||2 < (Mx,x) < e ||x|p for all x G ^. (3.11) 
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For the proof of Lemma [H we need the following general facts included for the conve- 
nience of the reader in form of the next two lemmata. 

Lemma 2. Let {X, (-, ■)) be any inner product space with norm \\ ■ |p := (-, ■). Moreover, 
assume that x,v G X fulfil the conditions \ \v\ \ = \ \x\\ and (x, v) > 0. Then, for P = j^v, 
we have 




X — v\\ < 11(7 — P)a;|| < ||a; — w 



Proof: First, we note that it suffices to show the equivalence for \\v\\ = \\x\\ = 1 and 
V ^ X. The second inequality is clear due to the fact that P is an ^2-orthogonal projector. 
For the inverse estimate, let = {x,v) > 0. We then have ||x — f |p = + — 

2ao = 2 — 2aQ and ||(/ — P)x|p = 1 — Oq, which together gives 



||(/-P)X||2 


1 + ao 


\v — a; P 


2 



□ 



Lemma 3. Let (■,-)2 be two equivalent inner products on a Hilbert space Ti, i.e. 

they induce equivalent norms in Ti, G a symmetric operator with respect to (■,-)2 (i-^- 
{Gx,y)2 = {x,Gy)2 for all x,y G Ti.) and u G kerG. Then, if G is {■,■) 2- elliptic on the 
(■, ■) I- orthogonal complement of u, that is 

{Gx,x)2>c{x,x)2 (3.12) 

holds for all x with {x,u)i = 0, then G is also {■,■) 2- elliptic on the {■,■) 2- orthogonal 
complement of u, i.e. ^3.12i) holds for all x with {x,u)2 = with a possibly different 
constant c. 

Proof: Let x E TL fulfil {x,u)2 = 0. Decomposing x = Pix + x"*"^, where PiX = u 
denotes the ^-orthogonal projector, we obtain Gx = Gx"*"^ and therefore 



{Gx,x)2 = {x^\Gx)2 = {Gx^\x^')2 >\\x^'\r2 ~ \\x 



II- 



To prove the assertion, it remains to see that Hx-*"^!!! ~ ||a;"'"^||2 = ||3;||2; where x^ 



{u,u)2 

be a multiple of u such that ||w||i = and {x,v)i > 0, giving 



_ (^.")2 ^ _ ^_ This latter fact follows from Lemma [2] above, where we can choose v to 



||x-^l||l ~ llx-t^lli ^ \\X - V\\2 > - P2){X - V)\\2 = \\x-^'^\\2. 

This confirms the assertion. □ 
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Proof of lemma [T], For the ellipticity, we show the claim for all x with (x, u)a = 0; 
setting (■, ■)! = (■, ■)a in Lemma [3] proves then the coercivity on Vq. Indeed, for all such 

||x||2 

X, we have \\Ji\t > A, by assumption [21 and therefore 

((A-AC)x,x) = ||x||i-A||x||^ (3.13) 

li-A||x||^ (3.14) 

> :i-^7||x||^ + HA||x||2,-A||x|pc (3.15) 

(3.16) 
By 

< ((A-AC)w,w) < (Aw, w) < r||wf , (3.17) 
the boundedness of M on Vq follows. 

□ 



l|x||i-A||x| 


|2 

Ic 


A-A„ „2 


A„ 


A " 




A-A „ „2 
A 






A-A „ „2 




A 





Lemma 4. Let $:=(! — P)(I — aM.) : Vq ^ Vq denote the matrix in the left part of 
1^3. 8\) . Then setting the parameter a := 2((1 — A/A)7 + r)~^ in the Richardson method 
gives 

||*||yo-vb < < 1- 

Proof: Due to the continuity and ellipticity of M on Vq, we have 

(l-ae)||x||2 < ((I-aM)x,x) < (l-a^)||x||^ V x G K,. (3.18) 

The same inequality applies to $ owing to the fact that P is symmetric and therefore 

((I-P)(I-aM)x,x) = ((I-aM)x, (I-P)x) = ((I-aM)x,x), V x e K). 

We now let P = < 1. Then, choosing a := = 2((1 — A/A)7 + F)"^ (recalling the 
choice of Q and 6 from the proof of Lemma [1]), equation fl3.18p gives —(3 < x, x) < /3 
and therefore ||$||yo^vb < /5 < 1- D 



3.2.3 Estimates for Rayleigh quotients and residuals 

It remains to estimate the Rayleigh quotients in the right hand side of equation (13.81) . 

Lemma 5. Let x G i2{I) with ||x|| = 1 and as before let A(x) denote the corresponding 
Rayleigh quotient. Then we have 

9F 

A(x)-A < — A(x) ||5(x)^||2, (3.19) 

7 
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where the constants in inequality Ii3.19\) are those from the norm estimate Ii2. 2\] and 
^"""(x) = (I — P)(x — u) as above. 

Proof: Suppose that u^; is an eigenvector belonging to the smallest eigenvalue A of the 
eigenvalue problem (12.11) which is normalized in such a way that | lu^^l | = | |x| | and (x, u^;) > 
0. Straightforward computation yields 

A(x)-A = ^^(llx-u^lli-Allx-u^lPc). (3.20) 



c 



Keeping the normalization of x and (12. 2p in mind and since there are no absolute 
values involved there, we infer from (I3.20p that 



1 ,, ,,o l|x 



|2 



A(x)-A < ---T2-||X- U^IIa = -p— rT2-||x- U^.||a < ,, ..^ ||x - U^. 





X 


2 
A 


7l 


X 


|2 

Ic 



oil ^ llA II 112II ^ llA — II 112II 1 1 A 

il^llc ll^llc 7l|x|lc 

*X ^^C^ II I 1 2 A / \ I I I 1 2 

= II^-u^-IIa < -A(x)||x- u^-ll . 

7 7 

Since Lemma [2] says that 

||x-u,|| < V2||(I-P)x|| = v^||(I-P)(x-u,)||, 
the assertion (I3.19P follows as well. □ 

Remark 3. As before let 6 = <5(x) := x — u and fix < a < 7/2r. // ||(5(x)-'-|| < a/q, 
one has for the corresponding Rayleigh quotient 

A(x)<A(l-^)^'=:i^. (3.21) 

Taking e.g. a = 7/4r one has K = 2A. 

Proof: From Lemma[5]we infer that A(x) < A(l - ^||(5(x)^||2)"^ < A(l - ^)"^ □ 

An immediate consequence of (I3.2ip can be stated as follows. 

Remark 4. Under the property of Section\^ the operator C is coercive on any sufficiently 
small fixed angular neighborhood of the direction u, i.e. for K from Ii3.21\) one has 

(Cx,x)>7//s: u;/ienewr ||x|| = 1, ||(u-x)^|| < Va. (3.22) 

The next observation is that the orthogonal error components are controlled by the 
residuals. To this end, note that 

M±: Vo^ Vo, X (I - P)Mx 

is bounded and elliptic so that ||M7^|| is a bounded positive number. 



13 



Lemma 6. For any x G £2(1) with ||x|| = 1, let r(x) := (A — A(x)C)x be the corre- 
sponding residual with respect to the Rayleigh quotient A(x) := (ind, as before, let 
x"*- = (I — P)x so that in particular 6 (x.)^ = (I — P)(x — u) = (I — P)x. Assume that for 
the constant a from Remark\^ 

Then for M ■= 2\\Wl]}\\ one has 

2KVCc 
7 



||5(x)^|| < M||r(x)^|| < M||r(x)|| < m(||M|| + — -^||5(x)^||) ||5(x)^||. 

(3.24) 



Proof: Since 

r(x) = M(x - u) + (A - A(x))Cx = M(I - P)(x - u) + (A - A(x))Cx, 
we conclude, on account of Lemma [5] and Remark [3l that 

||r(x)|| < ||M|| \\5\x)\\ + '^-^^\\8\x)\\\ (3.25) 

7 

which confirms the upper estimate in f l3.24p . 

As for the lower estimate, recall that by the definition of M_l, 

Mi5(x)^ = (I - P)(A - A(x)C)(5(x)^ + (I - P)(A(x) - A)C(5(x)^, 

giving 

||Mji-i<5(x)^|| < ||M^(5(x)^|| < ||(I-P)(A-A(x)C)5(x)^|| + ||(A(x)-A)C(5(x)^||. 

(3.26) 

Now straightforward calculations yield 

(A - A(x)C)<5(x)^ = r(x) - (A - A(x)C)Px = r(x) - j^(A - A(x))Cu, 
which gives 

||(I-P)(A-A(x)C)5(x)^|| < ||r(x)^||+Cc|A-A(x)|. (3.27) 
Invoking now Lemma [5l f l3.19p together with Remark [3] and f l3.26p . yields 

||Mli|ri||5(x)^|| < ||r(x)^|| + ^^^(1 + ||5(x)^||)||5(x)^|p. (3.28) 
Hence, whenever 

^^^(1 + ||5(x)^||)||5(x)^|| < \m-^\\-\ (3.29) 
which is the case when (13.230 holds, we obtain 

||5(x)^||<2||Mli||||r(x)-^||, (3.30) 
which is lower estimate in (13.241) with M := 2||M2"'^||. □ 
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3.2.4 Convergence of the scheme 



In summary, the above observations allow us to establish the following local convergence 
properties of the Richardson eigenvalue iteration. 

Theorem 1. Let the damping parameter a in the basic algorithm MINIT from Section 
\3.1\ he chosen according to Lemma^ so that the hound (3 for ||$||yQ_»yg in Lemma ^ 
satisfies jS < 1. Furthermore assume that the initial error satisfies 

\\S^\\<mm\a'/',^}=:d, (3.31) 



2dK 

where c := 2arCc/7- Then the following statements hold: 
a) For ^ := (/3 + l)/2 < 1 one has 

llCill <eil'5;^l|, (3.32) 

i.e. one has monotone linear error reduction with respect to the i2-norm. Moreover, for 
any e > there is an such that for n > one can choose ^ < (3 + e in liS. 32\] . i.e. 
one has an asymptotic error reduction of a rate < (3. The orthogonal error components 
8^ also converge to zero in the A-norm. 

h) There exists a uniform constant C such that the error of the Rayleigh quotients A(x,„) 
is hounded hy 

A(x„)-A < C\\8i\\^ ^ 0. (3.33) 
c) There exist constants M, M' > and uq & N such that 

\\Si\\ < M\\ri\\ < M||r„|| < M'\\S^\l n > no- (3.34) 
Moreover, there exists C < 1 o,nd ni E N such that 

||rn+i|| < Clknll, n>ni. (3.35) 
Proof: Recalling (13. 5p . (13. 6p and employing Lemmata H] and [5] provides 

ll^n+lll < ll^llvo^K) ll<^nll + ^ K^n)\\Sn\\^ 

< (/? + cA(x„)||(5;^||)||(5;^||=:e||5;[||, (3.36) 

with c := 2aTCc/'y- It remains to show that ^ < 1 for \\Sq\\ sufficiently small. In fact, we 
know from Remark [3] that 

A(xo) < K, (3.37) 

provided that 6o = <5(xo) satisfies \ \Sq\ \ < y/a (cf (13.311) ). Moreover, by (I3.36P for n = 0, 
we have 

115^11 < {P + cK\\di\\)\\di\\. 
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Now note that (13.311) implies that ci^"||5^|| < (1 — /9)/2 so that 

^ + cK\\5^\\ < ^^=--i < 1, (3.38) 

which yields 

\\5^\\ < A(xi) < J^. 
One easily shows now with the aid of Lemma and (I3.36P inductively that 

\\SUi\\<U\Si\l in-=P + cK\\Sn\\<i<^. A(x„)<i^, neN. (3.39) 

This guarantees monotone convergence of the iteration. It is obvious from the above 
arguments that with 1 1<5^| | — >■ the values of ^„ drop down to (3, showing the first assertion. 

Convergence in the A-norm is now an immediate consequence of (12.21) and (13.321) 
provided that 5q satisfies (I3.3ip . This proves a). 

The convergence of the Rayleigh quotients A(x„) now follows immediately from (13.191) 
and (I3.39p . giving (13.331) . confirming b). 

Concerning c), we have seen above that under the assumption (I3.3ip the estimates 
(I3.39P hold. Thus, using the rough estimate a < 1 one obtains 

^^(I + II'^^^IDII^^II < ^l|M-||-\ n>no, (3.40) 
7 ^ 

where 

no : 



'log(7/(8rKCc||M 



-1 1 



log^ 

Hence, for n > uq, the hypotheses of Lemma [6] are satisfied for x = x„, which provides 
the lower estimate of (I3.34p again with M = 2||M2^|| as in Lemma O The upper estimate 
is also an immediate consequence of Lemma E] and (13.391) . 

As for the remaining claim, we find that 

||l"n+l|| = ||(A — A„,+iC)x„+i|| 

= ||x„+i||-i ||(A-A("+i)C)(x„-arO|| 
= ||x„+i|ri ||(A-A(")C)x„ 

-a{A - AC)r„ + (A^") - A("+^))Cx„ - a{X - A("+^))Cr„| |. 

On account of Lemma O (I3.19p . (I3.2ip . and the fact that | |x„_|_i| | > 1 for the unnormalized 
iterates, we obtain further 

||r„+i|| < ||(I-«(A-AC))|| ||r„|| + ^^!^||5;[|p(2||x„|| + «||r„||) 



7 

< [P + Y (2 + a||r„||j||r„| 
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where M is from fl3.24p and where we have also used Lemma HI the continuity of C and the 
estimate (13.241) just proven before. Recall that we have /3 < 1, thus showing the assertion 
for ||r„|| small enough. The fact that ||r„|| indeed becomes small follows from the upper 
estimate in (13.241) and (I3.32p as well. □ 



4 Adaptive strategies for the eigenvalue problem 

There are two points that we would like to stress from the start. First, all the above 
considerations are independent of the index set 1 underlying the space ^2(2^) being finite or 
infinite as long as our assumptions on A, C are satisfied. If the basic Richardson iteration 
is considered on a fixed finite dimensional space, i.e. if we consider a fixed discretization of 
the original problem (11.21) like (11.71) . the proposed scheme may not be the most favorable 
one since it offers at best linear error reduction. Nevertheless, the property (12. 2p reflects 
the fact that one has been able to precondition the problem sufficiently well in the sense 
that a fixed error reduction is achieved independently of the (fixed) size of X. 

On the other hand, if one is interested in solving the original problem (11.21) within a 
desired target accuracy e at possibly low cost (which means to understand how the cost 
depends on e when e tends to zero), the game changes completely. Following [CDD2j one 
may think of performing the iteration on the infinite eigenvalue problem ( 12. ip which is 
still equivalent to i \1.2\\ . In fact, solving (12. ip within target accuracy e provides a solution 
to (II. 2p . due to the norm equivalence (12.20 and the mapping property (II. 5p . that has up 
to a constant factor the same accuracy in the (continuous) energy norm. Now, of course, 
the matrices A, C are infinite so that even when the current iterate x„ has finite support, 
the next iterate x„+i cannot be computed exactly. The idea is therefore to compute each 
update only approximately within a suitable dynamically varying accuracy tolerance in 
such a way that, on one hand, the convergence to the exact solution is preserved, while 
on the other hand, the computational cost of each perturbed iteration is possibly low. 

In the following sections we shall carry out this program based on the above Richardson 
type iteration. It will be seen that, regardless of the order of the iteration, when applied 
to a fixed finite dimensional trial space, such an adaptive scheme may perform at an 
asymptotically optimal complexity. 

We will proceed in two steps. First, each iteration step requires the approximate cal- 
culation of the residual (A — A(x„)C)x„ which can be broken down to two types of tasks, 
namely: 

(I) to approximate matrix/vector products Ax, Cx (where now A, C are in principle 
infinite). 
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(II) to approximate the Rayleigh quotient A(x) = 

Of course, (II) can be reduced to (I) but there are interesting alternatives. We distinguish 
the two cases at this point because, depending on the situation, the accuracy requirements 
may be somewhat different due to the fact that relative accuracy tolerances are needed 
in (II) unless C has particularly favorable properties that essentially permit the exact 
calculation of (Cx, x). We shall return to this issue later in more detail. 

Let us therefore suppose at this point that we have for B G {A, C} a routine with the 
following property at hand: 

APPLY{B,x,r]) w such that 

I |Bx — w| I < rj. 

We shall postpone the actual description of APPLY to a later section. Given such a 
routine we shall first analyze which tolerances t] are needed to ensure convergence of a 
correspondingly perturbed iteration. As a second step we shall then discuss the complexity 
of a corresponding numerical realization. 



4.1 A perturbed iteration 

We shall now collect the main ingredients of a perturbed version of MINIT from Section 
13.11 assuming that the routine APPLY from above is available. The first one concerns 
the approximation of a given x G ^2(2^) by one with possibly short support. 

APPROX{:K,r]) z: produces for a given (finitely supported) x and any t] > a finitely 
supported z such that 

||x-z||<7/, #suppz < argmin||^_^ll<^/2#suppw. (4.1) 

This routine can be realized by replacing as many entries of the input as possible by 
zero, as long as the sum of their squares does not exceed 77^. One could of course use the 
tolerance 77 in the argmin which would mean to compute a best A-term approximation of 
X. This in turn would require exact sorting of the coefficients introducing an additional 
log-factor of the support size of x, see |CDDlj for details. Being content with quasi-sorting 
based on binary binning one can avoid the log-factor at the expense of a slightly larger 
support, see e.g. [B]. 

In addition to approximate matrix/vector products we need approximate Rayleigh 
quotients and hence scalar products. Recall that we have to perform such routines for 
normalized inputs ||x|| = 1, which will be henceforth assumed. A straightforward way to 
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approximate the scalar products would be as follows: 



SCAL(B,:si,T]) s: Given any rj > 0, the routine outputs a scalar s such that 

|(Bx,x) - s\<r] (4.2) 

as follows: 



• APPi?OX(x, r//(2 1 |B 1 1)) z; 
. APPLr(B,x,r//2)|,,pp, -.w; 

• s = (z, w) . 

Here APPLy(B, x, ?7/2)|guppz means that the output of APPLY can be restricted 
to suppz. To which extent this can be used to economize on the computational effort 
depends on the realization of APPLY and will therefore be postponed to Section 14.31 

It is easily seen that s given as above does satisfy (14. 2p . Note that 1 1 A| | := sup||x||=i 1 1 Ax| | < 

r. 

We also remark that one can think of different ways of approximating (Bx, x). For 
instance, one could decompose x by determining Xj of smallest support so that | |x — Xj 1 1 < 
2^7^1/2^ say, noting that xj = for J := \\ log2 77^/^|]. Then, setting zj := Xj — Xj+i 

J 

J](z„ w,), := APPLY{B, x, e,)|,,pp,^,, , j = . . . , J, (4.3) 

i=o 

is an approximation to (x, Bx) of the order 77, provided that X]j=o ^''^i ~ ^7^^^- The appar- 
ent advantage is that highly accurate matrix applications need to be computed only on 
typically small supports of coarse Zj's. It will be shown in Section 1473) that such a strategy 
does offer asymptotic savings. Of course, as long as the matrix/vector products Ax, Cx 
have to be approximated within a similar accuracy tolerance anyway in the iteration, one 
might as well stick with the simpler version of SCAL shown above which we will do for 
the time being. 

Given the routine SCAL we can proceed to approximating the Rayleigh quotients 
A(x). We shall devise a routine: 



RAYL{'x,ri) — > A: Given rj > and any finitely supported :k with ||x|| = 1, RAYL outputs 
A auch that 

|A(x) - A| < r/. (4.4) 
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In many cases it is even possible to apply C exactly at acceptable cost, e.g. when C 
is diagonal, see e.g. Remark [H In this case we need no further assumption on C beyond 
boundedness in ^2(^) and can simply take 

fl^rL(x,,) = ^^^%^4^. (4.5) 

(Cx, x) 

If, on the other hand, C can only be applied appoximately with the aid of some routine 
APPLY, the routine SCAL is slightly more involved. Denoting by cc '■= l/K the local 
coercivity constant for C from Remark [3l the following realization works. 

Lemma 7. Assume that ||(u — x)-*"!! < ^/a holds. Carrying out the following steps: 

• SCAL{A,^,r]Cc/2) ^ §a, 5CAL(C, x, c2,r//(6r)) ^ 5^; 

• A := Sa/sc, 



verifies (i^); provided that 

r] < min {7/cc, 31/ cc}- (4.6) 

Proof: To see the validity of (14.41) . let sa '■= (Ax, x), sc '■= (Cx, x) and note first 
that 

sa sa I , ~ X sa , ~ X 
— = — [sa-sa) —[sc-sc), 

Sc Sc Sc ScSc 

so that 

|A(x) - A| < p + J^f . (4.7) 

2sc scsc er 

Now note that under the assumption fl4.6p one has (recalling that ||x|| = 1) 

\sa-Sa\ < — = < ^ < -SA, 

SO that Sa < < Likewise, upon using (I3.22p . one obtains for rj satisfying (14.61) that 
l-sc — sc| < sc/2 so that sc > sc/2 > cc/2. Combining these estimates yields 

Sa ^ er 



scSc 2c^ ' 

which, on account of (14.71) . confirms the claim (14. 4p . □ 



Remark 5. On account of Remarkl^ the tolerance in the evaluation of the scalar product 
remains proportional to the target accuracy rj in the routine RAYL for either version 
provided that ||(u — x)-*"!! < ^/a. 
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We are now able to describe a routine for approximating the scaled residual a(A — 
A(x)C)x: 

RES{:K,ri) — > r^: Given any rj > and any x with ||x|| = 1, the routine RES outputs a 
finitely supported vector such that 

\\rr, - a{A - A(x)C)x| | < r/. (4.8) 

The routine RES can be realized as follows: 

1) RAYL{^,r]/{4Cca)) A; 

2) RES{x,r]) = a{APPLY{A,^,r]/{2a)) - XAPPLY{C,^,7]/{AXa))) 

Of course, when the approximate evaluation of the scalar product involving A simply 
relies on the application of APPLY one only needs to carry out this latter routine once 
in RES with respect to the minimal tolerance required in RAYL and in step 2) of the 
above scheme. 

Remark 6. From now on we shall always assume that the parameter a in the basic 
Richardson iteration is chosen according to Lemma^ ensuring that (3 <1. 

The following lemma will help to understand the perturbed Richarson iteration. 

Lemma 8. Assume that the approximation u to u satisfies \\6{vl)-^\\ < e < d, where d 
is the constant from liS. 31\) . and let ^ be given by Ii3. 38\) . Then, setting rjj := (1 — ^)e2~^ , 
Vo := u, and doing for j = 0, 1, 2, 3, . . . 



one has 



and 



\j - RES{\j, r]j) Vj+i, 

||5(v,)^|| <^", J = 0,1,2,... 
^{^r^)^\\<ee//3, j = 0,1,2,... 



(4.9) 

(4.10) 
(4.11) 



Proof: Setting r^ := RES{\j,ri), we have, by definition, 

Vj+i-u = - u - ar(vj) + (ar(vj) - r^J 

= Vj — u — a(A(vj — u) — AC(vj — u)) + a;(A(vj) — A)Cvj 
+ (ar(vj) -r,,J 

= (I - aM)(v, - u) + «(A(v,) - A)Cv, + (ar(v,) - r,J. (4.12) 
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Therefore we obtain, on account of Lemma HI Lemma 0, and the accuracy of r^^. 
||5(v,+0^ll < /3||5(v,)^|| + ^^S:A(v,)||5(v,)^ir + r/, 

= (/5+^^y^Mv,)ll'5(v,)^ll)ll^(v,)^ll+^r (4.13) 

Now recall from Remark [3] that A(vj) < K a,s long as | |<5(vj)-'-| | < ^/a < d (which is, in 
particular, valid by assumption at the initialization step vq = u). 

From (13.391) we then know that 

||(5(vi)^|| < (/5 + ^i^i^||5(u)^||)||5(u)^|| + (l-0^-<ell'5(vo)^|| + (l-0^-<^-^ 

where ^ < 1 is given by (13.381) . Hence, we conclude that | |<5(vi)-'-| | < e. One easily 
concludes inductively that (14.101) holds (actually with strict inequality for j > 1). 

More precisely, repeating the reasoning in (I4.13p . we obtain upon elementary calcula- 
tions and using that ^ = (1 + /3)/2, 

ii5(v,+o^ii < j^e-v. < (4.14) 

which was to be shown. □ 



Now recall from Lemma [HI that, once the orthogonal error component drops below 
another possibly smaller threshold 



2IV FA' C M7^ 



1/2 



(4.15) 



the orthogonal error component of the approximate eigenvectors behaves essentially as 
the residual. 



Remark 7. If 5 < d at most 



mo 



log 



I3S 



log^ 



iterations of the form suffice to provide an approximation x^^ to u satisfying 

||5(x„,J^|| <min{ci,5}. (4.16) 

We are now ready to formulate the main adaptive eigensolver. On account of Remark 
[7] we shall assume for simplicity without loss of generality that the initial guess already 
satisfies the somewhat more stringent accuracy tolerance (I4.16p . 
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MINIEIG(A,C,eH (A(£),u(£)) 

(i) Initialization: Choose Eq < nim{d,d} (sec p.3ip . (|4.15p ) and ||xo|| = 1 s.t. j|(5(xo)^j| < Eq; 
set e := Eq, u xg; 

(ii) Iteration block: set u — > v, (1 — ^ r], j = 0; 

do 

RESiv,!-]) ^ r^; v - r,, ^ v; v/\\v\\ ^ v; t]/2 ^ rj, j + I ^ j; 
until j >\og{ciP)/ log C or 7/+ ||r,,|| <cie/M for a fixed a < min {2/(5eo), (V3 + 5/2)-i2-3/2}; 
(Hi) Coarsening: do w = APPRO X{-v,3ci£/y/2) ; 

normalize u = j^', if < £, stop and set A(e) = RES{vL,e), u(£) := u; 
else e/2 e; go to (ii) . 



Note that a block of perturbed iterations (ii) is interrupted by a coarsening step (iii) 
as soon as a threshold criterion is met. This criterion involves actually two alternatives 
that will both be seen to be met after a uniformly bounded finite number of steps in (ii). 
One of the stopping tests is an a-posteriori test based on the numerical residual and the 
rationale is that it may actually be met earlier than the other test which is based on the 
bounds in (14. lip which might be too pessimistic, in particular at later stages when the 
reduction constants decrease. The role of the coarsening step is to control the complexity 
of the overall scheme in a similar way as in adaptive schemes for operator equations (see 
e.g. [CDD211CDD3]). 

Theorem 2. The scheme MINIEIG(A,C, e) terminates for any given e > after 
finitely many steps and outputs an approximate eigenpair {X{e),u{e)) with u normalized 
satisfying 

||(u-u(£))^|| <£, |A(5)-A| <£, (4.17) 
where (A, u) is an exact ground state eigenpair. 

Remark 8. Given u{e) satisfying the first estimate in lji4-i'A ) '^^ can approximate A within 
a tolerance proportional to instead of ((^. j7| j hy applying RAY L{\i{e) e^) which of course 
requires a correspondingly higher cost due to the the required more accurate applications 
of SCAL and ultimately of APPLY , see Section 

Proof of Theorem [2t To analyze the effect of the various perturbations of the exact 
iteration suppose that := u is the output of the coarsening step after the fc-th cycle 
through the coarsening step (iii), i.e. at this stage we have e = 2~^eo- Moreover, let Vj be 
the result after j perturbed iterations in the iteration block (ii) starting with vq = u = u^. 
The corresponding tolerance rj = rjj is then given by r/ = (1 — C,)b2~^ = (1 — ^)eQ2~'^~^ . 
Invoking Lemma [HI ensures that ||5(vj)"'"|| < e and | |<5(vj)-'"| | < e^-' / [3 and hence that the 
error reduces at least by a fixed rate. Hence, the stopping criterion in (ii) is met after at 
most J := [log(ci/5)/ log^] steps giving 

||5(vj)^|| <ci£. (4.18) 
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Moreover, the initialization ensures that Lemma E] apphes to the iterates Vj, which means 
that ||5(vj)-'-|| is controlled by the residuals. Thus 

||5(v,)^|| < M ||r(v,)|| < Af(||r,J| + r/,)- (4-19) 

On the other hand, also by Lemma [6l 

||r,,|| < ||r(v,)||+r/, < M'||5(v,)^|| + r/, < M'^eV/^ + r/,. 

Thus \ \r,y\ \ +rij will also drop below the threshold Cis/M after finitely many steps, and in 
fact, possibly earlier than the alternative step bound J, if the bounds in fl4.1ip are overly 
pessimistic. 

In summary, the input v in (iii) at this stage satisfies ||(u — v)-*"!! < cie. Assuming 
without loss of generality that (u, x) > 0, we can invoke Lemma [2] to conclude that 

||u° - v|| < V2cie =: a, (4.20) 

where u° := u/||u||. Since by (iii), we have w = APPi^OX (v, 3(t/2), we obtain by 
triangle inequality ||u° — w|| < 5(t/2. Moreover, by definition of the routine APPROX , 
we have (w, v — w) = so that ||w|| > (1 — {3a/2)'^Y^'^. Thus, the normalized vector 
u ■= w/||w|| satisfies with (1 - (3(t/2)2)-V2 =-. I + g 

||(u — u)"*"!! < I |u° — u| I < I |u° — w| I + I |w — u| I 

So- 11 11 / 1 ,\ 5(7 
< h w - — --1 < h w 



2 V||w|| J - 2 V(l - (3cr/2)2)i/2 

< — + g. 

- 2 ^ 

Noting that for 6 < 1/2 one has (1 - 62)-i/2 < i + 2b/ V^, we conclude that g < y/Sa, so 
that 

||(u-u)^|| < (V3 + 5/2)cr = v^Ci(V3 + 5/2)e < £/2, (4.21) 

by our choice of the constant ci. Thus, in summary we have shown that after a uniformly 
bounded finite number of perturbed iterations in (ii) with initial accuracy e, one branches 
into (iii) whose output is either sufficiently accurate or serves as input for (ii) with im- 
proved accuracy 6/2. Hence the algorithm terminates after finitely many cycles through 
(ii), (iii), namely as soon as 2~''eo < e. □ 



4.2 Complexity estimates 

It remains to analyze the computational complexity of the above scheme MINIEIG. The 
subsequent analysis follows similar lines as used before in connection with adaptive solvers 
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for operator equations. We formulate an ideal benchmark which describes the minimal 
cost needed to achieve a desired accuracy tolerance e for an approximate normalized 
eigenvector. This lower bound is simply the number of entries needed in any finitely 
supported sequence to approximate u within accuracy e (or equivalently the normalized 
u° within a fixed factor). This naturally leads to the notion of best A^-term approximation 
that we briefly recall first. 

Let Sfc denote the set of all x G ^2(2^) which have at most k nonvanishing entries. 
Then 

(j/v(x) := inf I |x — z| I 

is the error of best A^-term approximation in ^2(2^)- Then 

|x|_4s := sup A^*crAr(x) 
AfeN 

is a ( quasi- )seminorm and 

:= {x e £2(2^) : ||x|U^ := ||x|| + |x|_4. < 00} 

is a (quasi-) Banach space. Thus, for s > 0, the unit ball of is the set of all those 
sequences in ^2(2^) whose error of best A^-term approximation decays at least as A^~** 
and hence a compact set in £2(1) ■ Another way to view this is the following: In order to 
approximate a given x G with accuracy e it takes at most A^^ = £:~^/*|x|^/ entries to 
do so, and in the worst case over the whole unit ball of A^ the necessary order of entries 
is exactly e~^^^. In what follows this relation 

accuracy e <-> degrees of freedom e""^''''! ■ |_4s 

reflecting s-sparsity of elements in -^2(2^) will be a central orientation in the subsequent 
developments. 

Let us pause to mention that A"^ can also be characterized as a weak ip space. In fact, 
denote by x* the nonincreasing rearrangement of x G ^2(2^)? i-e. x*_^i = \xj^_^_j^\ < x* = \xj.\, 
j = 1, 2, . . .. Then ip{I) is comprized of all those x G ^2(2^) for which 

|x|c(x) := supn^/^x* < 00. 

n>l 

It is not hard to show that (see |De Vorej ) 

A"" = i = s + - 

Moreover, it is easy to see that ip C ip but for any q < p one has ip G ig, so that s-sparse 
sequences are almost just ip summable sequences with s and p related as above. 

The first key ingredient of the analysis is the following coarsening lemma that explains, 
in particular, the role of step (iii) in MINIEIG, see [Ul ICDD3j . 
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Lemma 9. Assume that v G for for some s > and suppose that x G i2{I) is any 
finitely supported sequence in £2(1) such that ||v — x|| < s. Moreover fix b > and set 

w := argmin||^_^||<(i+^^#supp z. (4.22) 

Then there exists a constant C depending only on b and s when s — > such that 

||w|U. < C||v|U., (4.23) 

and 

1 1 

#supp w < C||v||_^,£: ||v - w|| < (2 + &)£ < ||vm.(#supp w) ^ (4.24) 
hold uniformly in e > 0. 

In other words, thresholding a given finitely supported approximation at a slightly 
higher tolerance than the accuracy of approximation provides essentially a best TV-term 
approximation to the (possibly unknown) approximand. 

As mentioned before, strictly speaking the cost of determining w is essentially 
(#suppx)) log(#suppx)). However, at the expense of a slightly worse target accuracy 
than (1 + b)e one can get away with quasi-sorting based on binary binning at a computa- 
tional cost that stays proportional to (#suppx)). For simplicity, we shall therefore ignore 
the log-factor in what follows and use that APPROX can be realized in linear complexity 
of the input size. 

Remark 9. Let us denote again by Uk the result of the k-th application of step (Hi) in 
MINIEIG. Then, if the ground state u belongs to A'^ for some s > we have 

llufclU^ < ||u|U., #suppufc < (£o2-^)-^/l|u||i//, ||u°-Ufe|| <£o2-', (4.25) 
uniformly in k eN. 

Proof: The last relation in (14.251) has been already established in the proof of Theorem 
[21 see ( I4.2ip . The rest is then an immediate consequence of Lemma [91 □ 



Hence to estimate the computational complexity of MINIEIG{A,C,e) depending 
on e as e ^ 0, it remains to bound the computational work in each block (ii). 

Proposition 1. Assume that for some s > one has u G A'^ , and that r^^ := RES{\,'r]) 
satisfies 

W^vWa^ < c'(||u|Us + ||v|Us), 



#suppr„ #fiops(r,) < Cr/"!/^ (| |u| 1^// + | |v| 1^//) 



(4.26) 
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holds for some constant C independent of t]. Here #flops(r,;y) denotes the number of 
arithmetic operations needed to compute r^. Then the output {X{e), u{s)) of MINIEIG{A, C, e) 
satisfies in addition to ( [y^.i7| ) 

#suppu(£), #flops(u(£)) < Ce-'/'\\u\\](:, \\u{e)\Us < C||u|U., (4.27) 

for some constant C independent of u and e. 

Proof: Applying if necessary the coarsening lemma to the initial guess (starting from a 
correspondingly slightly higher initial accuracy) we have I xoUs < |u|_4s, ^supp xq < ^o^^"*- 
Suppose now that at the k-th. stage of (ii) the input u = \ik-i satisfies for some constant 
C independent of k 

#suppu, #fiops(u) < Ce~^''{\\vL\\'i! + ||u||^//). (4.28) 

By assumption 04.261) . the same estimates hold for the intermediate iterates \j in (ii) 
(with Vq = Ufc_i) with constants however now depending on j. As shown in the proof of 
Theorem [2] each block (ii) has at most a uniformly bounded number J of iterations, so 
that fl4.28p still holds for the input to step (iii) with some constant C = C{J), J being 
the upper bound of the number of iterations in (ii). As pointed out in Remark [9], the ap- 
plication of APPROX produces then u = Ufc satisfying again f l4.28p for Sk = Sk-i/'^ with 
a constant coming from the coarsening lemma that is independent of k. Thus summing 
the computational cost of each block (ii) and using a straightforward geometric series 
argument confirms the claim. □ 

Thus, it remains to verify f l4.26p for the approximate residuals. According to the ingre- 
dients of the approximate residuals provided by the routine RES, given in the previous 
section, the key issue here is the efficient approximate application of the matrices A and C 
through the routine APPLY. It is well known by now that for a wide range of (local and 
global) operators and suitable wavelet bases the corresponding operator representations 
in wavelet coordinates are nearly sparse. A precise formulation of this property that is 
fulfilled in many concrete cases, reads as follows, see e.g. |CDD4l [DHS] . 

s*-compressibility: Lets* be a positive real. B^* denotes the set of matrices (overXxX) 
with the following properties: B G B^* if for every /c G N, there is a matrix B^, with at 
most 1^ entries in each row and column satisfying ||B — Bfc|| < C2"'^* ak, where (afc)^Q 
is an £i-sequence of positive numbers. Elements in B^* are called s*-compressible. 

One can show that any B G B^* maps the ^''-spaces boundedly into themselves 
for s < s*, cf. |CDD1] . Section 3. For elements in B G B^* a concrete realization of 
APPLY{G, - , ■) has been developed in [CDDl] whose properties are given in the following 
lemma. 
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Lemma 10. Assume that B G B^* . Given a tolerance 6 > and a vector x with finite 
support, the algorithm AP P LY (B , :si, 6) produces a vector w which satisfies 

||Bx- w|| < 6. 

Moreover, for s < s* one has: 

(i) The output vector w satisfies 

l|w|U» < ||x|U»; supp w < ||x||^,5~i 

(ii) The number of entries ofB to be computed to obtain w is < ||x||j^s(5^~;. 
For the proof again see jCDDlj . 

If for s < s* the number of arithmetic operations needed to compute w does not exceed 

1 _ 1 

C| |x| \^sS » + 2 supp X (using quasi-sorting instead of exact sorting which would entail an 
additional log-factor), the matrix B is called s*- computable, [S| IGSj . For the verification 
of s*-computibility for a wide class of operators, see [Si IGSj . 

s*-sparsity: The matrix B is called s*-sparse if there exists a scheme APPLY satisfying 
the properties listed in Lemma [73 whose computational complexity remains proportional 
to the output size. 

Clearly s*-computable matrices are examples of s*-sparse matrices. It is important to 
note though that there are further important examples. Recalling the form of C in remark 
m the matrix A may actually be the product of several matrices. An APPLY-scheme for 
such products can easily be obtained by composing individual APPLY -schemes designed 
e.g. for compressible matrices, see [CDD2] for the treatment of least squares formulations. 

Another important case concerns the application of C in a fairly general setting. 

Actually C represent the inverse C = of an (elliptic) operator S in the sense of 
(11. 5p . Thus the APPLY-scheme for C may just mean the adaptive approximate solution 
of an operator equation for an 7i-elliptic operator S : Ti ^ Ti.'. In many cases its output 
has been shown to satisfy the properties in Lemma [10] and thus gives rise to an s*-sparse 
APPLY-scheme. 

In summary, it is important to keep in mind that the actual realizations of APPLY 
for A and C may be completely different but should satisfy the properties of Lemma [TOl 

Property 3. The matrices A and C are s* -sparse for some s* > (For C this is trivially 
the case when C can be applied exactly in linear time). 
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The main result of this paper can now be formulated as follows. 

Theorem 3. Assume that A and C are s* -sparse for some s* > and that the parameters 
a, (3 are chosen according to Lemma^ Then for any s > 0, the scheme MINIEIG{A, C, e) 
after finitely many steps produces an approximate eigenpair (A(e),u(e)) with normalized 
u satisfying 

\\{u-u{e))^\\ <e, |A-A(£)|<e, (4.29) 

where (A, u) is the exact ground state solution of Ii2.1\) . 

Moreover, if u E for some s < s* , then one has 

#flopsu(£), #suppu(£) < s-'/'\\u\\'j!, \\u{e)\Us < \\u\Us, (4.30) 

where the constants are independent of s and u but depend only on s when s approaches 
s*. 

Proof: f l4.29p has been already shown in Theorem [21 To prove the rest of the claim we 
employ Proposition [1] which requires confirming the property (14.261) . Towards this end, 
recall from Section |4?T] that the realization of RES{'x,ri) requires the approximate appli- 
cation of A, C within tolerances that are uniformly bounded from below and above by 
fixed multiples of rj (see (14.81) ) as well as the computation of RAY Lijs., rj / (Aa)) . Now, this 
latter routine requires the evaluation of the routine SCAL for A and, unless C can be 
applied exactly, for C. By Remark [5] and Lemma [71 the tolerances needed in the SCAL 
routines remain uniformly proportional to the target accuracy in RAYL{:k, rj/ (4a)) which 
is proportional to t]. Since all these routines rely on the scheme APPLY with respect to 
tolerances proportional to rj, the relations (14.261) follow from Lemma [TOl This finishes the 
proof. □ 



4.3 Possible quantitative improvements 

The scheme analyzed above should be viewed as one possible realization of an adaptive 
strategy. To achieve quantitative improvements of schemes of the above type one might 
try to exploit the fact that the (exact) Rayleigh quotient exhibits essentially the square 
accuracy of the corresponding approximate eigendirection. We shall only sketch such 
strategies whose details would essentially follow analogous lines as discussed above. 

Recall that the evaluation of the scalar products requires the tightest accuracy tol- 
erances (although still proportional to the target accuracy of RES). Thus there are two 
angles that might help to reduce computational complexity, namely: 

a) Trying to speed up the calculation of scalar products; 

b) Reducing the number of calls of RAYL. 
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4.3.1 Fast evaluation of scalar products 

As for a), the naive approach outlined above rests on the accurate calculation of a com- 
plete matrix/vector product although one computes at the end only a single number. 
We have already mentioned a possible approach in (14.31) that might help reducing the 
computational complexity. 

Recalling that an approximate eigendirection of accuracy e would give rise to a (pre- 
cise) Rayleigh quotient that approximates A within accuracy of the order of (see Lemma 
|5]), a first natural idea is to postprocess the result of MINIEIG{A, C, e), satisfying (14.291) . 
so as to obtain an approximation to A of order e^. In fact, in view of Lemma O it would 
make sense to compute 

X*{e) = RAYL{u{e),e^), so that |A*(e) - A| < (1 + C)^^, with C := 21/^/7, 

see (I3.19P and (13.2ip . However, using the simple version of the routine SCAL described 
below (14. 2p . the computational complexity could be of the order of e"^/'' when u G A^. 
Let us now point out that this cost can be reduced significantly by employing more refined 
versions of SCAL along the lines of (14. 3p . To this end, suppose that B is an s*-compressible 
matrix and recall e.g. from [CDDlirPHSj the following key idea of constructing an APPLY 
scheme for an approximate matrix/vector computation. Fix any s < s* . Given C > and 
V of finite support, let v<^ := APPROX{y, Q and set 

V-i(C) := V - V(^, Vj(C) := V2is-<. - V2(i+i).-<^, j = . . . , J(C) := [logad |v| |/C)/s] . 

Then, one has 

J2 v,(C) = V, ||v_i(C)|| < C, l|v,-(C)ll < (1 + r)2n, J = 1, • • • , AO- (4.31) 

i=-i 

Moreover, it is shown in [DHSj that (when Bj are the compressed versions of B G B^* 
from the definition of s*-compressibility) 

:= ^BjVj(C) satisfies ||Bv - w^-|| < C'C, (4.32) 

j=0 

for some uniform constant C independent of (. For simplicity we shall work with C = 1 
which can always be arranged through the definition of the Bj or by replacing ( by c( in 
the decomposition of v. 

The announced improved version of SCAL is based on the following observation. 
Proposition 2. Assume that B G B'^* and s < s* is fixed. Setting 

00 

:= aj{l + 2')-^2-'^^/^ , j = 0, . . . , J(v^), e_i := a_iv^, ^ aj = 1, (4.33) 
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where the summahle (to one) coefficients aj are again from the definition of s* -compressibility 
and have algebraic decay, and defining 

s(5,v):= 5^(v,(v^),w,^), (4.34) 



where the w^^, are given by ^4-32 ), one has 

|s(5,v) - (Bv,v)| < 5. (4.35) 

Moreover, whenever supp v < Cr^/^'^lvl^// for some s < s < s* , then the number 
of operations ops{v,6) needed to compute s(5, v) is bounded by a constant multiple of 

s+2s 

^|v|^s5"2i j , where the constant depends only on C . 
Proof: By (14.311) we can write 

\s{6,^r) - (Bv,v)| = | (v,(v^),w,^ _ Bv)| < ^ a,S < S, 

j=-i j=-i 

where we have used Cauchy-Schwarz and (14.331) in the last step, confirming (I4.35p . 

As for the work count, we shall estimate now first the number of operations needed 
to compute a single entry of w^. Now note that (since v G for any s > 0) for Aj : = 
suppVj(C) one has 

#A, < (2^-^"C)-^/lv|^//. (4.36) 

Further, recall that each row of Bj has at most 2^ entries, so that the computation of one 
entry of a contribution BjVj((^) takes, in view of (14.361) . at most 2^ operations as long as 
j < j* when j* = j*{C) is the largest integer for which 



(S + S)-M0g2 (- 



c 



(4.37) 



Thus the computation of a single entry of the partial sum Ylj=o'^j^jiC) takes the order 
of 

2i* < (4.38) 
Likewise the number of computations required for a single entry of the remaining sum 



Ej=riO B,v, (C) is, by (gSS]), of the order of 

J 

c'^'\^\](:J2^~''^' < c'^'\y\](!2'^*'/' < (4.39) 



j=j* 

Therefore 



#ops (for computing one entry of w^) < ( ^Z'-'^"'"'^-' |v|^i (4.40) 
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We shall now estimate the work required by the computation of (v(j), w^^.). Note first that 

#suppv_i(v^) < r^/^lvl^//, #suppv,(v^) < 2-^"^>ri/2^|v|i//. 

Hence, using f l4.40p with ( = ej, the computation of {\^j{^/6),w^.) in (14. 341) takes, in view 
of (14.341) . the order of 



#suppv,(v^)|v|^/i^+')eTi/(^+^") < a/+'{l+2') — \Y\p''2-^'"/'^''+'"^6 '^\^+'J (4.41) 



operations. Summing over j and recalling that the aj decay polynomially completes the 
proof. □ 



Corollary 1. Assume that the hypotheses of Theorem\^ are valid. Given 5 > 0, let u(e) 
he the output of MINIEIG{A,C,e) and let \*{e) = RAYL{u{e),e'^) where RAY L is 
based on a version of SCAL derived in an obvious manner from ^.34\) - Then one has 

\\-\*{s)\<{l + 2VKh))e\ (4.42) 

and, whenever vi G A'^ for some s < s, the computational complexity i^ops{\*{e)) of \*{e) 
remains bounded by 



i^ops{\*{e)) < \n\p-'h =V=+V, (4.43) 
where the constant is independent of u and s. 

Proof: The estimate (14.421) is an immediate consequence of (13.191) and the first rela- 
tion in (I4.29p . The complexity estimate, in turn, follows from Proposition [2] applied to 
V = u{e) with 5 := together with fOOD . □ 

Since g{s) := (s + 2s)/(s + s) increases in s and g{s) = 3/2, the computational 
complexity of computing X*{e), and hence the smallest eigenvalue A within a tolerance of 
order e"^, grows at most like which is of course much better than the cost e"^/* that 

would result from applying the original simple version of SCAL. In fact, when s is very 
small one almost recovers cost of e~^^^ needed to approximate u G within tolerance e. 

Note also that the coarsened versions of v needed in the computation of w^^ are 
essentially the same as those in (I4.34p . so that they can be reused. Nevertheless, this 
version of RAYL is quantitatively still more involved as the original simple one based 
on approximating Bv at the respective target accuracy. Since these vectors are needed 
anyway in the course of MINIEIG it seems preferable to use the latter more elaborate 
version only at the end once u(£:) has been obtained. 
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4.3.2 Modified iterations 



As for b), another option is to modify the ideal iteration itself, again trying to exploit 
the fact that the approximation rate of the (exact) Rayleigh quotients is faster than that 
of the eigendirections. We only give a very rough sketch of the idea. Instead of the final 
accuracy one could run MINIEIG first with some intermediate accuracy e outputting 
again \i{e) satisfying (14.291) . One could then fix A = \*{e), set xq := u(e) and iterate 

Xj+i = Xj - a(A - AC)xi, 

so that, by straightforward calculations, one obtains 

(x,+i - u)^ = (x, - u)^ - $(x, - u)^ + a{\ - A)(I - P)Cx„ 

and hence by recursion 

||(x.-u)^||<^..+ "+f:f')-^ |C||, (4.44) 

Thus, in principle, after |log£:| steps one has quadratic accuracy of (xj — u)-*- without 
further applications of RAYL. However, the iterates Xj are no longer orthogonal to the 
residuals := (A — AC)xj. Observing that (fi,Xj) = (A(xj) — A)(Cxj,Xj), we have 



|x,+i|p = ||xi||2W||ri||2-2a(A(x,)-A)||xi||^ > \\^i\\^W\\ri\\^-2a{\{^,)-\)\\C\\\\^i\\^, 



so that we cannot guarantee any more that the iterates have increasing norms. Neverthe- 
less, since |A(xj) — A| is expected to be of the order and ||xo|| = 1, there is still a fixed 
positive constant 6 > 1/2 say, so that after J ~ | loge| steps, one still has ||xj|| > 6, so 
that renormalizing the Jth iterate preserves quadratic accuracy. This way, one expects to 
catch up quadratic accuracy of the approximate eigendirections without any intermediate 
computation of Rayleigh quotients. One can then continue such a block of iterations with 
initial guess xj/||xj|| that approximates u now within a tolerance of the order of so 
that a corresponding approximate eigenvalue can be computed within tolerance with 
the aid of the above fast computation of scalar products, etc. We leave the details to the 
reader. 



Concerning optimal line search (13. ID and subspace acceleration techniques, we mention 
that this requires the solution of a two-dimensional (or more generally small) generalized 
eigenvalue problem. However, setting up the corresponding matrices requires the compu- 
tation of scalar products. Since, generally, we cannot compute these quantities exactly, a 
careful assessment of the perturbation effects would be necessary. Also, the possibilities 
for realizing enhanced accuracy at reduced cost mentioned above might be relevant in this 
context. 
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5 Remarks concerning an application - 
The Schrodinger equation 



We conclude with some brief comments concerning a scenario that has actually motivated 
part of this work. 



We consider the time- independent electronic Schrodinger equation Hip = Eqi/j, where 
H : H\{R^ X {±|})^) H~\{R^ x {±|})^) is the Hamilton operator of the molec- 
ular TV-electron system under consideration, and Eq denotes the lowest eigenvalue of H 
corresponding to the ground state of the system. The spectrum of H is real and bounded 
from below by a certain /i G M (cf. e.g. |Yser] and references therein), so that the shifted 
Hamiltonian H' := H — fiS satisfies the conditions (II. 4p and (II. 5p of Section [TJ Letting 
$ = {(puW ^ ^} be an orthonormal basis of L2((M^ x {±i})^), we have to solve the 
eigenvalue equation 

He = EqC (5.1) 
where the entries of the "matrix Hamiltonian" in this formulation are given by 

In view of our discussion in Section [1], the problem is posed in ^2(2^)) which gives the 
so-called complete CI formulation for the Schrodinger equation. To obtain a Riesz ba- 
sis for , we can utilize a reference operator JF = ^^^^ -Fi, where the Fi : H^{E? x 
{±i}) H~^(^ X {±|}) are single particle operators with a complete eigenvalue ba- 



sis. The self consistent Fock operator from the Hartree Fock equations, see e.g. Helg 
can be modified to fit the present purpose. This operator is also bounded as an operator 
H^{{m? X {±\Y) i/~i((M3 X {±\}Y) and its spectrum is bounded from below, so 
that JF' := T — iiE satisfies (11.51) . In this case the eigenf unctions and the eigenvalues 
Oy of JF, T' can easily be computed from those of the single particle operator Fj, 

If we choose 

^ = {i^y := a'y'^4>y\v ^T) 

this basis satisfies (II. 8p . The generalized eigenvalue problem ( 12. ip then resembles the 
symmetry-transformed variant of the Schrodinger equation with 

C := ((a;^0„a70,,)),yex and A := C^H'C^ = C^(H-/iI)C^ (5.2) 
where the conditions (12. 2p and (12.30 are valid. 
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With the above notation, ||x||c then gives the L2((M^ x {±i})^)-error of the wave 
function := J^^ei^'^'^i' ~ Tliv&^i^'t'y For estimating the H^{{^ x {±|})^)-error of 
the wave function note that the canonical H^{{M? x {±i})^)-norm is equivalent to 
the norms induced by the inner products ((., .)) := {H'-, and (JF'-, defined by the 
Hamiltonian and the reference operator respectively, which gives 

||x||i = (C^H'C^x,x) = (H'C^x,C^x) 

and 

||x||2 = ||C5x||c-i = ||c||c-i = llv^ll^' - llV'll/fi- 

This serves also as a typical example where C is not coercive on £2{^) (because ||.||c is 
equivalent to the L2-norm), cf. Remark [H The norms ||.||a and ||.||2 — ||-||a on ^2(2^) can 
now be used to estimate the convergence of the CI solution with respect to the iJ^-norm. 
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